{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Satellite Ground Contacts\n",
"\n",
"This tutorial computes communication contact intervals between a satellite and a set of ground stations over a single day.\n",
"\n",
"A satellite is \"in view\" of a ground station when it rises above a minimum elevation angle (here, 5 degrees). For each second of the day, we compute the satellite's position relative to each ground station in the North-East-Down (NED) frame, determine the elevation angle, and identify contiguous time intervals where the satellite is visible.\n",
"\n",
"The example uses the Landsat-7 satellite and three ground stations: Svalbard, Alice Springs, and Sioux Falls."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scienceplots # noqa: F401\n",
"plt.style.use([\"science\", \"no-latex\", \"../satkit.mplstyle\"])\n",
"%config InlineBackend.figure_formats = ['svg']\n",
"from datetime import datetime, timedelta\n",
"\n",
"# The TLE for landsat-7\n",
"tle_lines = [\n",
" \"0 LANDSAT-7\",\n",
" \"1 25682U 99020A 24099.90566066 .00000551 00000-0 12253-3 0 9992\",\n",
" \"2 25682 97.8952 129.9471 0001421 108.5441 14.5268 14.60548156329087\",\n",
"]\n",
"landsat7 = sk.TLE.from_lines(tle_lines)\n",
"\n",
"# The minimum elevation for a contact\n",
"min_elevation_deg = 5\n",
"\n",
"\n",
"# The date to compute ground contacts: April 9, 2024\n",
"date = datetime(2024, 4, 9)\n",
"# Any array of times representing every second of the day\n",
"time_array = np.array([date + timedelta(seconds=x) for x in range(86400)])\n",
"\n",
"# Get satellite positions in TEME frame (pseudo-inertial) via SGP4\n",
"pTEME, _vTEME = sk.sgp4(landsat7, time_array)\n",
"\n",
"# Get ITRF coordinates (Earth-Fixed) by rotating the position in the TEME frame\n",
"# to ITRF frame using the frametransform module\n",
"pITRF = np.array(\n",
" [q * x for q, x in zip(sk.frametransform.qteme2itrf(time_array), pTEME)]\n",
")\n",
"\n",
"# Setup some ground stations\n",
"ground_stations = [\n",
" {\"name\": \"Svalbard\", \"lat\": 78.2232, \"lon\": 15.6267, \"alt\": 0},\n",
" {\"name\": \"Alice Springs\", \"lat\": -23.6980, \"lon\": 133.8807, \"alt\": 0},\n",
" {\"name\": \"Sioux Falls\", \"lat\": 43.5446, \"lon\": -96.7311, \"alt\": 0},\n",
"]\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compute Satellite Positions\n",
"\n",
"Propagate the TLE with SGP4 to get TEME positions at every second of the day, then rotate to the ITRF (Earth-fixed) frame. Define the ground stations by their geodetic coordinates."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Contact Detection Algorithm\n",
"\n",
"For each ground station, the algorithm:\n",
"\n",
"1. Computes the satellite position in the station's **North-East-Down (NED)** frame by subtracting the station position and rotating with `qned2itrf`\n",
"2. Computes **elevation** as the arcsine of the upward (-Down) component of the normalized NED vector\n",
"3. Finds all times where elevation exceeds the minimum threshold\n",
"4. Groups consecutive qualifying indices into individual contact intervals\n",
"5. For each contact, records the start/end times, duration, range, elevation, and heading"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Contact Summary\n",
"\n",
"Display all contacts for the day sorted by start time. Higher maximum elevation generally corresponds to better signal quality."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Contact Detail\n",
"\n",
"Plot the range, elevation, and heading over time for an individual contact pass."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def calc_contacts(ground_station, pITRF, time_array):\n",
" \"\"\"\n",
" Compute contact times for a single ground station given satellite position in Earth-fixed frame\n",
" \"\"\"\n",
"\n",
" # Create an \"itrfcoord\" object for the ground station\n",
" coord = sk.itrfcoord(\n",
" latitude_deg=ground_station[\"lat\"],\n",
" longitude_deg=ground_station[\"lon\"],\n",
" altitude_m=ground_station[\"alt\"],\n",
" )\n",
"\n",
" # Get the North-East-Down coordinates of the satellite relative to the ground station\n",
" # at all times by taking the difference between the satellite position and the ground\n",
" # coordinates, then rotating to the \"North-East-Down\" frame relative to the ground station\n",
" pNED = np.array([coord.qned2itrf.conj * (x - coord.vector) for x in pITRF])\n",
"\n",
" # Normalize the NED coordinates\n",
" pNED_hat = pNED / np.linalg.norm(pNED, axis=1)[:, None]\n",
"\n",
" # Find the elevation from the ground station at all times\n",
" # This is the arcsine of the \"up\" portion of the NED-hat vector\n",
" elevation_deg = np.degrees(np.arcsin(-pNED_hat[:, 2]))\n",
"\n",
" # We can see ground station when elevation is greater than min_elevation_deg\n",
" inview_idx = np.argwhere(elevation_deg > min_elevation_deg).flatten().astype(int)\n",
"\n",
" # split indices into groups of consecutive indices\n",
" # This indicates contiguous contacts\n",
" inview_idx = np.split(inview_idx, np.where(np.diff(inview_idx) != 1)[0] + 1)\n",
"\n",
" def get_single_contacts(inview_idx):\n",
" for cidx in inview_idx:\n",
" # cidx are indices to the time array for this contact\n",
"\n",
" # the North-East-Down position of the satellite relative to\n",
" # ground station over the single contact\n",
" cpNED = pNED[cidx, :]\n",
"\n",
" # Compute the range in meters\n",
" range = np.linalg.norm(cpNED, axis=1)\n",
"\n",
" # elevation in degrees over the contact\n",
" contact_elevation_deg = elevation_deg[cidx]\n",
"\n",
" # Heading clockwise from North is arctangent of east/north'\n",
" heading_deg = np.degrees(np.arctan2(cpNED[:, 1], cpNED[:, 0]))\n",
"\n",
" # Yield a dictionary describing the results\n",
" yield {\n",
" \"groundstation\": ground_station[\"name\"],\n",
" \"timearr\": time_array[cidx],\n",
" \"range_km\": range * 1.0e-3,\n",
" \"elevation_deg\": contact_elevation_deg,\n",
" \"heading_deg\": heading_deg,\n",
" \"start\": time_array[cidx[0]],\n",
" \"end\": time_array[cidx[-1]],\n",
" \"max_elevation_deg\": np.max(contact_elevation_deg),\n",
" \"duration\": time_array[cidx[-1]] - time_array[cidx[0]],\n",
" }\n",
"\n",
" return list(get_single_contacts(inview_idx))\n",
"\n",
"\n",
"# Calculate all the contacts\n",
"contacts = [calc_contacts(g, pITRF, time_array) for g in ground_stations]\n",
"\n",
"# Flatten contacts into 1D list\n",
"contacts = [item for sublist in contacts for item in sublist]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Convert to pandas dataframe for nice table display\n",
"\n",
"import pandas as pd\n",
"\n",
"data = pd.DataFrame(contacts)\n",
"data.sort_values(by=\"start\", inplace=True)\n",
"data.reset_index(drop=True, inplace=True)\n",
"# Get nicer column names for display\n",
"data.rename(\n",
" columns={\n",
" \"max_elevation_deg\": \"Max Elevation (deg)\",\n",
" \"duration\": \"Duration (s)\",\n",
" \"start\": \"Start (UTC)\",\n",
" \"end\": \"End (UTC)\",\n",
" \"groundstation\": \"Ground Station\",\n",
" },\n",
" inplace=True,\n",
")\n",
"data.style.hide(\n",
" subset=[\"timearr\", \"range_km\", \"elevation_deg\", \"heading_deg\"], axis=1\n",
").format(\n",
" {\n",
" \"Max Elevation (deg)\": \"{:.1f}\",\n",
" \"Start (UTC)\": lambda x: x.strftime(\"%H:%M:%S\"),\n",
" \"End (UTC)\": lambda x: x.strftime(\"%H:%M:%S\"),\n",
" \"Duration (s)\": lambda x: x.seconds,\n",
" }\n",
")\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Plot one of the contacts\n",
"contact = data.iloc[5]\n",
"\n",
"fig, axes = plt.subplots(3, 1, figsize=(10, 7), sharex=True)\n",
"\n",
"axes[0].plot(contact[\"timearr\"], contact[\"range_km\"])\n",
"axes[0].set_ylabel(\"Range (km)\")\n",
"axes[0].set_title(f\"Landsat 7 to {contact['Ground Station']} on {contact['Start (UTC)']}\")\n",
"\n",
"axes[1].plot(contact[\"timearr\"], contact[\"elevation_deg\"])\n",
"axes[1].set_ylabel(\"Elevation (deg)\")\n",
"\n",
"axes[2].plot(contact[\"timearr\"], contact[\"heading_deg\"])\n",
"axes[2].set_ylabel(\"Heading (deg)\")\n",
"axes[2].set_xlabel(\"Time\")\n",
"\n",
"fig.autofmt_xdate()\n",
"plt.tight_layout()\n",
"plt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.14.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}